Quantum Statistics and Thermodynamics in the Harmonic Approximation 



O 
(N 

o 

Q 

oo 

Oh; 



> 
m 



X 



J. R. Armstrong, N. T. Zinner, D. V. Fedorov, and A. S. Jensen 
Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark 

(Dated: December 30, 2011) 

We describe a method to compute thermodynamic quantities in the harmonic approximation 
for identical bosons and fermions in an external confining field. We use the canonical partition 
function where only energies and their degeneracies enter. The number of states of given energy 
and symmetry is found by separating the center of mass motion, and counting the remaining states 
of given symmetry and excitation energy of the relative motion. The oscillator frequencies that 
enter the harmonic Hamiltonian can be derived from realistic model parameters and the method 
corresponds to an effective interaction approach based on harmonic interactions. To demonstrate 
the method, we apply it to systems in two dimensions. Numerical calculations are compared to a 
brute force method that is considerably more computationally intensive. 
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I. INTRODUCTION 

The many-body problem cannot be solved exactly for 
realistic interactions and systems. Numerous approxima- 
tions have been formulated and applied over the years. 
One of the problems is that many particles require a large 
Hilbert space to allow for the many possible correlations. 
In fact, the space typically grows exponentially with the 
particle number. The necessary reduction of the Hilbert 
space to obtain a tractable problem has to be accompa- 
nied by a correponding transformation of the interaction 
which in turn has to be used in the smaller space. This 
is the method of effective interactions in reduced spaces. 

A different approach starts by an approximation ap- 
plied to the interaction itself. In many physical fields, the 
two-body interaction is a complicated entity, and there- 
fore it is desirable to reduce this complexity. One such 
direction is to perform a harmonic approximation of the 
interaction Hamiltonian which leaves all terms (one- and 
two-body) as a polynomial of at most second order in 
the coordinates of all particles. This yields an exactly 
solvable model for the many-body system [l|, 0] . How- 
ever, to utilize such an approach, the input parameters 
must be adjusted to reproduce properties of realistic sys- 
tems for low particle numbers. This could be two-body 
binding and scattering information and structural expec- 
tation values. 

This sort of approach has been pursued in nuclear 
and condensed-matter physics for many years [1, 0] as 
a source of exact insight into many-body problems [1] 
when numerics are either intractable or hard to inter- 
pret. More recently, cold atomic gases have emerged as 
an arena and testing ground for various models due to 
the large degree of controllability in experiments . In 
particular, interaction can be controlled through inter- 
atomic Feshbach resonances tTj. However, these systems 
are almost exclusively studied in the presence of an ex- 
ternal confining potential which can often be modelled 
by a (possibly deformed) harmonic oscillator. For exper- 
imental setups with very deformed traps, the geometry 
can in fact be manipulated to study effective one- or two- 



dimensional dynamics. Very recently, it has even become 
possible to explore the few-body limit with trapping of 
just a few atoms near the degenerate regime [1]. 

At the moment, the physics of ultracold atomic gases 
is, however, dominated by short-range interacting neutral 
atoms Q. In addition, since the samples are often ex- 
tremely dilute, the details of the potential at short-range 
does not matter, and only low-energy scattering informa- 
tion is important Q. The replacement of the potential 
by a regularized zero-range potential is therefore often 
justified. In the presence of a harmonic confinement, the 
solution of the two-body problem with a zero-range po- 
tential was presented by Busch et al. [l^. This solution 



was subsequently confirmed in experiments jll| . More 
recently, this model has been used as a starting point for 
addressing the few- and many-body problem in a con- 
trolled manner in both nuclear and cold atom physics 

The zero-range approximation scheme is one approach 
to effective interactions. Here we pursue a different one 
via harmonic expansion of the interaction Hamiltonian. 
System that can be described in some regimes of parame- 
ter space by harmonic Hamiltonians are extremely useful 
due to their solvability. Of course, as mentioned above, 
the parameters have to be meningfully related to realistic 
physical systems. This can be done through fits to two- 
body properties as described previously in Ref. . Here 
we take a parametric approach and study the details of 
the harmonic approximation for identical quantum parti- 
cles as function of the two relevant frequencies for such a 
system; the (internal) interaction frequency and the (ex- 
ternal) trapping frequencies, both of which are described 
by (isotropic) harmonic terms. Deformation can be easily 
implemented and will not be considered here. 

The objective of the present paper is to describe a 
method for calculating thermodynamical quantities for 
systems of identical bosons and fermions, a problem of 
relevance for any subfield of physics concerned with quan- 
tum mechanical behavior of multiparticle system in equi- 
librium. This requires full access to the partition function 
and the thermodynamic quantities that can be derived 
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from it. The quantum statistics of ensemble of particles 
is needed here, and we present a general method for ob- 
taining such information for harmonic Hamiltonians that 
works for any number of particles in principle. 

There is great interest in the universal thermodynam- 
ics of strongly interacting Fermi gases [isj . and mea- 
surements have determined the equation of state of the 
system [3]. Those studies are mostly concerned with 
trap-averaged quantities that are less sensitive to sud- 
den changes in thermodynamic behavior as expected near 
phase transitions. However, recently it has become pos- 
sible to determine local properties of the gas by density 
fluctuation measurements [15|. Quantities that involve 
derivatives of the energy or pressure like heat capacity 
and compressibility can thus be used to study for instance 
superfluid properties of both fermionic [Tcj and bosonic 
systems [17j . In addition, the same experiments are also 
able to explore effects of changing the dimensionality of 
the system through optical lattice potentials. 

In the current presentation we focus on extending th 
harmonic approximation into the statictical regime, and 
more particularly on the technicalities of the method used 
to achieve this. To demonstrate our method we consider 
the case of a two-dimensional system of identical bosons 
or fermions for different ratios of interaction to exter- 
nal trapping frequencies. The paper is organized as fol- 
lows. In SeclUwe briefly recap the harmonic method and 
then discuss the partition function in the harmonic ap- 
proximation. The most important part of the method is 
the counting of the degeneracies in the partition function 
which we discuss in extended detail. The thermodynamic 
quantities are briefly outlined and a discussion of how to 
related the parameters in the harmonic Hamiltonian to 
the model by Busch et al. is given. Practical numer- 
ical results are obtained and discussed in Sec. HIT] for a 
system of cold atoms confined to two spatial dimensions. 
Both bosons and fermions are explicitely investigated. 
Finally, give a summary and an outlook in Sec. IIVI 



II. METHOD 

We consider a system of N identical quantum particles, 
fermions or bosons, that have two-body interactions, all 
of which are confined by an external field. All interac- 
tions are substituted by carefully chosen harmonic os- 
cillators. The solutions of the Hamiltonian can then be 
found explicitly and treatment of thermodynamic prop- 
erties is possible. We first present the necessary ingredi- 
ents of the model, and continue to develope the method 
for calculating the partition function. This includes the 
key discussion of how to obtain the degeneracies of the 
many-body states. 



A. Formulation of the model 

The iV-body Hamiltonian contains kinetic energy, one- 
and two-body interaction terms, all of which we assume 
can be written in terms of second order polynomials 
in the particle coordinates when proper choice of effec- 
tive interaction parameters have been made. The result- 
ing harmonic oscillator structure can be separated into 
terms related to each of the Cartesian coordinates and 
the Hamiltonian in one, two, and three dimensions is a 
sum of terms from each coordinate. For example, in the 
x-direction, the effective Hamiltonian is 

^ (f 1 ^ 

= -^Erf^+2"^E-o4 (1) 

fc=i « 

where m is the mass of the identical particles and /i ~ 
to/2 is the reduced mass. The one-body (external) har- 
monic trap frequency is wq, while the two-body interac- 
tion frequency between all pairs is denoted Wint (inter- 
nal). VxQ is a constant energy shift which is included for 
generality. Double counting of interactions account for 
the factor 1/4. For two and three dimensions, the total 
Hamiltonian \s H = + Hy and H = + Hy + ■ We 
will assume that the external harmonic trap is isotropic, 
extension to the deformed case are straightforward and 
will be addressed in future work. All particles are identi- 
cal and interact in the same way which means that Wint 
is independent of i and k. Also the one-body frequency, 
u!o, is independent of k. 

The Schrodinger equation for the Hamiltonian in 
Eq. ([T]) can be solved analytically for any particle num- 
ber (details can be found in [2] ) The pertinent properties 
in the present context are that the energies, Ej, of the 
excited states have oscillator character. Here Ej denotes 
the jth many-body excited state, i.e. the energy of all N 
particles when they are distributed in a certain way in 
the oscillator levels that come out of the diagonalization 
of the full Hamiltonian. In this paper we will concentrate 
mostly on the case of a two-dimensional system, and we 
therefore specialize to this case now. The energy of a 
many-body state, Ej, can be written 

Ej = hujr{N,ci + N-1) + hujo{Ncm + 1) , (2) 

k=l 

iVcm = «r + "r ' (4) 

where and n!^™ are quanta in the x-direction related 
to the relative and center-of-mass excitations with the 
former carrying the index k running over all N particles 
(and similar terms for the j/-direction). The total number 
of relative and center-of-mass quanta in the many-body 
state j are denoted iVioi and A'cm respectively. Note that 
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n^, n^™, riy, and n™ must all be specified to completely 
characterize the many-body state of energy Ej. How- 
ever, more than one set of individual quantum numbers 
can lead to the same energy Ej and this determines the 
degeneracy of the state, which we denote gj and discuss 
at length below. Ej and gj are all that is needed to com- 
pute the partition function and thus the thermodynamic 
properties of the system. 

The constants, fvjJr{N — 1) and fno^ in the formula 
for Ej, are from the two-dimensional zero-point energy 
of intrinsic and center of mass parts, respectively. The 
one-body frequency, wq, correspondings to center of mass 
motion. The other frequency, uJr, is iV — 1 times degen- 
erate for identical particles. It is given by the relation 
i 

2 2-^,2 ,'r-\ 

^r=^i„ty+Wo- (5) 

If the two-body interaction goes to zero (wint — >■ 0), 
LOr Wo J a-nd the resulting eigenvalues and eigenfunc- 
tions converge to those of the external one-body trapping 
potential. The solutions are N non-interacting identical 
particles occupying N states in a harmonic oscillator. For 
interacting particles, we obtain instead — 1 identical 
frequencies with value given by Eq. ([5]). The extension 
of the above discuss to three dimensions involves also 
the z-direction and is straightforward, although the ex- 
tra quantum numbers makes the degeneracy larger and 
the counting of many-body states more involved. 

With the harmonic Hamiltonian, the exact spectrum 
of excited many-body states, Ej, are thus available for a 
given number of particles. For the thermodynamic cal- 
culations, the natural choice is to use the canonical en- 
semble applicable for a definite number of particles. The 
canonical partition function is 

Z = Y^g,eM~E,lT), (6) 

3 

where gj is the degeneracy of the jth many-body state 
and T is the temperature measured in energy units (the 
Boltzmann constant, ks — !)• Once Z is obtained, all 
thermodynamic quantities can be derived from it. The 
most involved question is how to obtain gj for a given 
energy of the many-body system. This will now be dis- 
cussed. 



B. Symmetry constraints on degeneracies 

The only quantities in the partition function are the 
energies, Ej, and their degeneracies, gj. The energies are 
appropriate multiples of the harmonic oscillator frequen- 
cies introduced above. The degeneracies depend strongly 
on which symmetry restriction is imposed. Problems 
arise when all, or a group of, particles cannot be distin- 
guished because they are identical and allow occupation 
of the same states. We then have to count the number of 



(anti-)symmetric states for the indistinguishable particles 
of a given excitation energy, Ej, and the corresponding 
degeneracies, gj, of each many-body state of energy Ej. 
We now describe our procedure to obtain gj. It is based 
on knowledge of the non-interacting system, a recursive 
procedure, and the fact that when the two-body inter- 
action vanishes, the degeneracies remain the same as we 
will now explain. 

We start with N identical particles without any two- 
body interaction (ojint 0), and all in the same ex- 
ternal trap with frequency loq. This Hamiltonian will 
have N degenerate solutions of frequency ujq per spa- 
tial dimension (so 2N for the two dimensional case). 
The many-body solutions are the products of N com- 
binations of the single-particle harmonic oscillator wave 
functions. For identical fermions, the Pauli principle, or 
equivalently the antisymmetrization, requires that each 
of the single-particle states at most is occupied by one 
particle. In addition, if this condition of occupancy by at 
most one particle for each state is fulfilled, then one and 
only one antisymmetric state is uniquely constructed as 
the Slater determinant of the N singly occupied single- 
particle states. 

Importantly, the energy can be characterized by the 
sum, A'tot — N-rci+Ncm, of single-particle oscillator quan- 
tum numbers for all particles in the state Ej, as seen in 
Eqs. dS]) and dH). Assume that we have computed the 
number, 5non(-^tot), of properly symmetrized states of 
non-interacting identical particles in the external poten- 
tial. The quantity gnon(-^tot) can be computed in a sim- 
ple manner by trail-and-error, i.e. taking all possible per- 
mutations and testing which are completely symmetric 
and which are completely antisymmetric. This is, how- 
ever, notoriously difficult since the computational effort 
is exponentially increasing with particle number. Also, it 
does not distinguish between center of mass motion and 
relative motion and this needs to be disentangled. This 
can be done by recalling the symmetry of the center of 
mass as we now discuss. 

The center of mass motion is a symmetric mode in all 
permutations of the particle coordinates, simply because 
the center of mass coordinate is the sum of the individual 
particle coordinates. This motion is always symmetric, 
independent of the number of quanta in this mode, Ncm- 
Therefore, the relative motion determines the symmetry 
of the total wave function. The excited states consist of 
a combination of two distinctly different types of excita- 
tions. They are characterized by the number of quanta 
in the center of mass motion, A'cm, and the number of 
quanta in the relative motion, iVioi- 

We denote the number of properly (anti)symmetric A^- 
body states by gp(N^ci), which is the number of states 
that have the proper symmetry in terms of the rela- 
tive coordinates. It can be computed from knowledge of 
the total number of non-interacting many-body states, 
5non(A^rci), for aU values of TVtot = N^ei + N^m- To do so, 
one must subtract the number of states of different center 
of mass quanta from the total number of non-interacting 
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many-body states in a recursive manner, i.e. 

5p(Mot) = 5non(^rcl) - ^ (iAr^„5p(A^rol " ^cm) , (7) 

Afcm = l 

where the degeneracy, rfjVcmj of the center of mass 
state depends on the spatial dimension of the oscilla- 
tor; diVc„ = 1 for one, dAr^,, = iVcm + 1 for two, and 
dwc^ = (-^cm + l)(^cm + 2)/2 for three dimensions, re- 
spectively. The formula comes from the fact that we 
know only the total number of quanta, A'tot, but not N^^i 
and A'cin individually. We need to recursively subtract 
the states that correspond to JVcm — 1,2,..., Ntot as 
the formula prescribes. The factor cIn^^ is due to the 
degeneracy of the center of mass motion itself for given 
number of quanta, A^cm- Note that 5p(0) is well-defined 
and corresponds to a state with no quanta of excitation, 
i.e. the ground state. It can be zero or non-zero de- 
pending on N and whether the particles are bosons or 
fermions. Eq. ([7]) determines the degeneracy of an iV- 
body state in a harmonic system and is one of the main 
results of this work. 

The degeneracy can now be calculated iteratively from 
Eq. ([7]) for any number of quanta, TVtot from a start- 
ing point based on a non-interacting system of identical 
particles in a trap. These states must then be cither 
completely symmetric or antisymmetric with respect to 
interchange of all coordinates of any pair of particles to 
obey either bosonic or fermionic statistics. Now we add 
the two-body interaction on top of the external one-body 
potential. The external frequency, wq, remains in the 
spectrum corresponding to the center of mass motion, 
and an iV — 1 times degenerate frequency, w^, appears 
corresponding to relative motion. The absolutely crucial 
point is that the degeneracies do not depend on the two- 
body interaction strength which will only influence the 
value of the degenerate frequency. The counting remains 
completely unchanged, because the number of states of 
a given symmetry is a discrete number which does not 
change with continuous variation of the interaction. 

Mathematically, this corresponds to a continuous map 
(the scheme for counting degeneracies) from a continu- 
ous interval (the interaction frequency, w^) to a discrete 
set (number of states, gj, with given energy). Such a 
map must necessarily be a constant map. The underly- 
ing point here is that our splitting into A^cm and iVroi is 
done in a way that moves the states that are related to 
the relative motion when is increased from zero to its 
full value, while the states that are related to center of 
mass motion remain constant in energy since ujq remains 
a solutions. 

We confirmed this simple counting method by an elab- 
orate and computationally very slow brute force proce- 
dure which is much worse than computing .gnon(A^tot) 
above. The (anti-)symmetrization of each of the com- 
puted wave functions is performed by permutations of all 
the coordinates followed by extraction of a basis with a 
size equal to the smallest number of linearly independent 



states. For sufficiently small N this can be computed 
in reasonable time and comparison can be made. In all 
cases considered below we found perfect agreement be- 
tween the numerics and the procedure outlined above. 
The procedure above is, however, far superior since it 
needs only counting of states with given energy, not the 
full wave functions. 



C. Thermodynamic quantities 

We can now use Eq. ([7]) to get the number of center 
of mass states and the number of relative states as func- 
tions of the sum of corresponding quanta. Along with 
the frequencies, this is all we need to calculate the parti- 
tion function in Eq. ([5]). This simplification occurs since 
the structure of the wave functions do not enter in the 
partition function. 

The center of mass mode of frequency loq separates in 
Eq. ^ because both the exponential function and the 
degeneracy factorizes, that is gj = gas{Nrei)dNam^ where 
the center of mass contribution is analytical and given 
below Eq. We have 



Z{N,T) 

Zcm{N,T) 
Zrel{N,T) 



= Z,„,{N,T) X Zrel{N,T), 

^ f cxpi-E,s/{2T)) Y 
\l-cxp{-Eo/T)) ' 



(8) 



gas{Nrel) ex.p{-Nrelfu^r/T) , 



where D is the dimension, Eq = hojQ, and Egs — 
Eq + {N — l)hajr + Vq- In practice, the infinite sum 
over relative quanta has to be limited by a cut-off value 
which is chosen sufficiently high to achieve the required 
accuracy. Certain quantities are more sensitive to the 
cut-off than others, such as those with more derivatives 
of the partition function (heat capacity, compressibility, 
etc.). For all the results presented below we have checked 
convergence by increasing the cut-off and identifying the 
value of T below which the relevant quantities remain un- 
changed. From the canonical partition function, Eqs. ^ 
or ([8]) , we get immediately the basic quantities of energy, 
E and free energy, F — E ~ TS, that is 



F 



-TlnZ 



dT 



dE 
df ' 



(9) 



where S is the entropy, and Cv is the heat capacity |18l . 

In thermodynamic formulations of macroscopic sys- 
tems, both temperature and volume is usually external 
parameters. Both the energy, E, and heat capacity, Cy, 
defined above, are in principle the derivatives for a fixed 
volume, V (and also fixed T and N). In our case with 
an externally confined A'^-body system the volume is not 
the fixed quantity, but rather it is the external trapping 
frequency, luq, that is fixed. The derivatives above are 
thus taken for fixed ujq. Usually, variation of the volume 
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allows information about pressure and compressibility. 
Definitions of these quantities and connection to their 
thermodynamic counterparts require a precise transla- 
tion between loq and volume. Here we employ the sim- 
plest choice and define the volume via the length param- 
eter related to the external trap, b, which is defined by 
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We now have V = sjjb , where S2 = and 



S3 = 47r/3 are the two and three dimensional angular sur- 
face areas. Derivatives of any quantity, Q, with respect 
to volume can now be written 



dQ dQ dujQ 2ujo dQ 
'dV ~ d^~dV ~ ^IWdu^' 

For the pressure, P, we find 

p_ dF _ 2^0 dF 

which, from the definition in Eq. gives 



P = Pc, 



P 



DV duio 
2Tu}a d\aZ^^i 
DV dujQ ■ 



(10) 



(11) 



(12) 
(13) 
(14) 



The derivatives are easily worked out since they only con- 
tain duJr/dujQ = UJo/uJr- 

Continuation to the second derivative with respect to 
the volume produce the isothermal bulk modulus, Bt, 
defined as 



R 1/ 1/^^ 2c.oap 



(15) 



which is the reciprocal of the isothermal compressibil- 
ity, kt- Again two terms arise related to center of mass 
and relative degrees of freedom. Several terms appear by 
carrying out the two derivatives. This kind of compress- 
ibility, that of a response of a system to an applied pres- 
sure, is different to the compressibility of a self-bound 
system such as a nucleus [20|. For self-bound systems, 
the compressibility refers to resistance to density fiuctu- 
ations. This can be written as the second derivate of the 
energy per particle with respect to the Fermi momentum 
for fermionic systems. 

Before we present results for degeneracies and various 
thermodynamic properties, we discuss how one can fix 
the two-body interaction term to capture the properties 
of the realistic system which have interactions that are 
not harmonic. In the case of identical bosons, the proce- 
dure was already presented in Ref. [2,] but we repeat the 
arguments here for completeness. 

We assume that we are in the situation that the real 
two-body interaction is of much shorter range than the 
external trap length b above. In this case, we can use the 
zero-range model of Busch et al. [l3| to obtain the solu- 
tions in the trap. The effective two-body interaction that 
we want to use is a harmonic oscillator, which in general 
contains two parts; an oscillator frequency, ojint, and a 




FIG. 1: The ratio between cOr from Eq. ((5]) and the exter- 
nal trap frequency, ujo, as function of the scattering length, 
a, in two dimensions for the Busch model [l^l with s-wave 
interactions for bosons (upper) and for the p-wave interaction 
identical (i.e. spin-polarized) fermions (lower). In both plots, 
the different curves correspond to different particle numbers. 
From the bottom to top, the curves are for 3, 4, 5, 6, 10, 
20, and 30 particles. The length parameter is the trapping 
length, I — b, introduces in the text. 



shift, Vq. To fix these parameters we employ the con- 
ditions that the two-body harmonic interaction should 
reproduce the correct two-body binding energy of the 
ground state in the Busch model, and also the spatial 
extend of the wave function. The latter condition is im- 
plemented by calculating the root-mean-square radius, 
(r^), in the exact solution and fixing the oscillator param- 
eter to reproduce value. The energy shift is suqsequently 
tuned so as to also reproduce the two-body binding en- 
ergy. 

Clearly, since the Busch model uses a zero-range inter- 
action, the only scale left is the external trap frequency, 
cjQ- Therefore what we obtain from this is the quantity 
Wint/wQ. We ignore the shift from now on since it merely 
provides an overall shift in the A'^-body harmonic Hamil- 
tonian which is not of interest in the current paper. The 
next step is to insert this value into Eq. ([5]) to obtain 
cjr/wo, which will now naturally depend on N. In the 
upper panel in Fig. [T] we show the result of this procedure 
for the case of bosons in two dimensions (more detail for 
bosons in both two and three dimensions can be found 
in Ref. 0)- The results are plotted as function of the 
two-dimensional scattering length a. Note that when a 
is small, the ground state becomes strongly bound, while 
for large a it goes to the non-interacting limit. This is 
clearly reflected in the behavior of w^/wo- 

The result of this procedure agrees rather well with 
predictions from the well-studied problem of bosons in- 
teracting through a zero-range interaction 21| in the 
strong and weak coupling limits 0, also in the case 
where higher-order interaction terms are included [2^ . 
Recently, a very similar procedure has been used to study 
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polar molecules in layered system [23| where excellent 
agreement with exact methods has been found both for 
isotropic [13] and anisotropic potentials [IE] which have 
reasonably large potential pockets. Incidentally, one- 
dimensional dipolar system also turn out to have such 
potential pockets [26l | and applying the harmonic approx- 
imation to these systems is an interesting direction for 
future research. 

For identical (spin-polarized) fermions a complication 
arises from the fact that the original Busch model ap- 
plies to particles interacting through an isotropic s-wave 
interaction only. The Pauli principle dictates that s-wave 
interactions are zero for identical fermions, and therefore 
the interactions must have p-wave (or higher odd par- 
tial wave) character. In two dimensions, the equivalent 
of the Busch model with p-waves can be solved and the 
spectrum turns out to be very similar to the s-wave case 
except for (unimportant) shifts 27]. However, a further 
complication arises since the corresponding ground state 
wave function is singular at the origin in a manner that 
does not allow it to be normalized [28[ . This can be fixed 
by a properly defined scalar product as discussed for the 
three-dimensional case in Ref. [1^. We will not pursue 
this approach here but instead we note that any higher 
order structural average of the type (r") with n an integer 
is still perfectly convergent. We thus consider a higher 
order average, {r*)/{r'^), in order to fix the oscillator pa- 
rameter for identical fermions. We are thus tacitly using 
that the p-wave spectrum is very similar to the s-wave 
one and assume that the wave function should be so as 
well. Of course, it must be kept in mind that this is not 
at odds with the Pauli principle since the two-body wave 
function in the p-wave channel is still zero at the origin. 
In any case, the Pauli principle is exactly enforced on the 
A^-body system as discussed above. 

The results for Wr/wo are shown in the lower panel of 
Fig. [T] The similarity to s-waves and bosons is quite 
clear, and the only difference seems to be slightly lower 
overall values for the fermionic case. Below we will illus- 
trate our method by calculating thermodynamic quanti- 
ties for both fermions and bosons. Here the only thing 
that matters is the ratio w^/wo- One can then reverse the 
process and use Fig. [T] to obtain a corresponding value of 
the scattering length and thus compare to a realistic sys- 
tem with trapped bosons or single-component fermions. 
Since this is mainly a discussion of method, we will not 
dwell on this any further. An important thing to note, 
however, is that since the ratios w^/wq are similar for 
bosons and fermions, the results below will truly iso- 
late the statistical behavior coming from the exchange 
requirements in the different cases. 



III. THERMODYNAMIC RESULTS 

We now proceed to demonstrate our methods by com- 
puting the density of states and the thermodynamic 
quantities that have been introduced above. They de- 
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FIG. 2: The number of completely symmetric (bosons) and 
antisymmetric (fermions) states as function of excitation en- 
ergy for different numbers of particles in a two-dimensional 
harmonically trapped system. The particle numbers are indi- 
cated on the plots. Notice that the horizontal axis for bosons 
is Ai'j.Yi^, while it is N^^'^N^^ for fermions. 



pend on particle mass, dimension, number of particles 
and their quantum statistics, temperature, one- and two- 
body frequencies. If we measure all energies in units of 
hujQ, all results for given particles in D dimensions, de- 
pend only on the two ratios, T/huQ and uJr/ujQ. This 
implies that any kind of interaction model used to define 
the effective harmonic hamiltonian only has to provide 
these two quantities as pointed out above. 



A. Degeneracies 

First we discuss the degeneracies themselves which 
are found by the procedure discussed above, essentially 
through the recursive formula Eq. ([7]) . While the number 
of non-interacting (anti-)symmetric states increases ex- 
ponentially as function of excitation energy, the numbers 
can, however, be tremendously reduced by the symmetry 
requirements. A simple example indicating this tendency 
is that for bosons all particles are allowed to occupy the 
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lowest level with zero quanta, Ntot = 0. In contrast, a 
given number of fermions has a minimum sum of quanta, 
Ngs, in the ground state, that is Ntot > Ngs ^ 0. We 
show the numbers of given symmetry as functions of ex- 
citation energy in Fig. [2j 

Since we are dealing with harmonic oscillators, a given 
energy corresponds to a given total number of quanta. 
The starting point for the ground state is naturally 
-^cm — and A^tot — -^rci, which is larger for fermions 
than for bosons due to the requirement that at most one 
fermion can occupy each state. The corresponding higher 
degeneracy implies therefore that the number of states 
of given excitation energy is larger for fermions than for 
bosons because the latter do not need to obey the Pauli 
Principle. For bosons, the lowest excited state consist- 
ing of one quantum added to zero quanta in the ground 
state can only be a center of mass excitation, since one 
quantum in the relative motion would be antisymmetric 
in one pair of coordinates (it would have to be a state 
of p- wave/parity-odd symmetry). The lowest completely 
symmetric excited state of the relative motion appear for 
A^ioi = 2. In conclusion, there is always a gap in the exci- 
tation spectrum for bosons corresponding to 2liujr where 
fermions only have a gap of huJr- This will be reflected 
in the temperature dependence below. 

For bosons (upper part of Fig. [5]) the same degeneracy 
is found for low energy (small total number of quanta) 
for all particle numbers. In the upper plot of Fig. [2] the 

1 /2 

horizontal axis is N^^^ since we have found that this is 
a good measure for the excitation energy in the system. 
For fermions, one can give a qualitative argument for 
the exponential behavior which we discuss below. The 
similarity of degeneracies at low energy for bosons can 
be understood by the direct counting procedure where a 
given number of particles, iV, has to be distributed to 
add up to the total number of quanta, iVj-ei- First, iV — 1 
particles are placed in the lowest oscillator level, and the 
one remaining particle then has to be placed in the level 
with total number of quanta equal to A'rci- Then we 
move the single particle one step down to A^rci — 1, and 
simultaneously compensating by moving one particle one 
step up from the lowest level. We continue with these 
combinations until we have iVroi particles in the second 
oscillator level and all others in the lowest level. Adding 
one particle and repeating the counting process we realize 
that their is a one-to-one correspondence between the 
configurations of N and -I- 1 particles. Going from one 
to the other is simply by removing or adding one particle 
in the lowest oscillator level with zero contribution to the 
total number of quanta. 

For N bosons, the deviation from this universal curve 
starts for N^d = N + 1. The reason is again found by fol- 
lowing the counting procedure. The configurations with 
all particles in the lowest two levels are only possible 
when iVrci < N. This implies that we find fewer states 
when iVrci > TV + 1 for iV than for TV + 1 particles. The 
curves break away from the universal curve for increasing 
N when N.^i = N + 1. 




kgT/(nco„) 



FIG. 3: Helmholtz free energies of a two-dimensional system, 
divided by particle number. A'', for several boson (upper) and 
fermion (lower) systems, identified by the number of particles 
and frequency ratio, [N,u)r/uJo). The ground state energy 
is subtracted, and the energy unit is hwQ for both energies 
and temperatures. The ground state energy, i?ga, has been 
substracted. 



The degeneracies for fermions have very different be- 
havior, as seen in the lower part of Fig. [5] Notice that 
the horizontal axis is different in the two plots. We no- 
tice a regime of linear dependence which has the same 
origin as the exponential square root dependence of ex- 
citation energy of the free Fermi gas level density [4j. 
For two dimensions, the particle number dependence is 
roughly N^^^ as reflected on the axis in the fermion plot. 
This holds for intemediate excitation energies and for 
relatively large particle numbers. A qualitative under- 
standing of the behavior can be obtained in a manner 
following Ref. The density of single-partile states 

at the Fermi level in a two-dimensional harmonic trap 
is roughly gp {2Ny/'^/hijJo, while the typical excita- 
tion energy in the system is E* = fiwo-^rei- The many- 
body level density in this situation is then proportional 



to exp 



This explains the choice 



of horizontal axis for fermions in Fig. [2j The suggested 
linear dependence is not very clear but the assumptions 
are not well fulflUed. The requirements is that the excita- 
tion energy has to be sufficiently large to allow statistcial 
treatment and sufficiently small not to exhaust particles 
at the bottom of the potential. This is not true for the rel- 
atively small particle numbers that we must necessarily 
work with to make the problem tractable by both brute 
force symmetrization and the counting scheme developed 
here. 



B. Energy and free energy 

We are now in a position to explore the thermodynam- 
ics and to compare bosonic and fermionic behavior in de- 
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FIG. 4: Energies in two dimensions, divided by NksT, for 
several boson (upper) and fermion (lower) systems, identified 
by the number of particles and frequency ratios, {N,iUr/ijJo). 
The ground state energy, _Ega, has been subtracted from each 
curve. 



tail. The free energies per particle for a two-dimensional 
system are shown in Fig. [3] as function of temperature 
for different interactions and particle numbers, both for 
bosons (upper panel) and fermions (lower panel). Note 
that the ground state energy has been subtracted as it 
is not of interest here. The curves start out very flat at 
zero showing that the low-temperature dependence is at 
least quadratic. The decrease with increasing temper- 
ature is actually rather similar to a quadratic behavior 
with T. The behavoir is modified when T exceeds huQ 
and the center of mass mode becomes active. Notice that 
the dependence on particle number is relatively weak. 
This is in contrast to the interaction dependence where 
the free energy varies substantially more for the smallest 
two-body interactions. In addition, we note that fermion 
free energies are always smaller than boson free energies. 
This demonstrates that the free energy is entropy driven, 
i.e. the entropy term grows faster than the energy with 
temperature. Again the larger degeneracy of fermionic 
systems is seen by the lower free energies. 

Fig. 0] displays the energy per particle for bosons and 




2 4 6 8 10 

kBT/(no)„) 

FIG. 5: Entropies in two dimensions, divided by Nks, for 
boson (upper) and fermion (lower) systems, identified by the 
number of particles and frequency ratios, {N,u)r/ijJo)- 



fermions. We measure energy in units of NksT to iso- 
late the high temperature behavior which is 2NkBT for 
two dimensionsal harmonic system by the equipartition 
theorem. Again we subtract the uninteresting constant 
ground state energy. The increase from zero at zero tem- 
perature is rather steep for both types of particles. All 
curves continue to increase with temperature but much 
slower after a few units of huQ. The largest energy per 
NksT is for the system with the least particles and the 
smallest frequency ratio. This can be understood from 
the fact that for large w^/wo, only the center of mass 
modes (twice degenerate in two dimensions) are active. 
However, since we divide by iV, the values become smaller 
for larger N, but the behavior remains the same (self- 
similar lines for Wr/t^o = 10 in Fig. 2]). This is simi- 
lar for both bosons and fermions. We also observe that 
the equipartition limit at high temperature is reached 
for substantially larger temperatures, a clear sign of the 
interaction effects. 

The low-temperature behavior is, however, different for 
bosons and fermions. The low-energy 2fiu}r gap in the bo- 
son spectra arises due to a ground state with all particles 
with zero oscillator quanta of excitation. Then there are 
no states with one quanta of excitation. For fermions, 
there are only gaps. Therefore the energy exhibits 
flat regions at relatively low temperature, but these fea- 
tures are observed first for fermions and later for bosons 
due to the difference in energy gap which is related to 
the activation energy. This demonstrates a clear signa- 
ture of quantum statistics in the harmonic model, but 
which should be expected in generic interacting systems 
with identical particles. 

Next, we show the entropy, S, in Fig. [S] We note 
the increase of available states per particle as function of 
temperature goes from zero to values around a few times 
the temperature in units of HijJq. For lower ratios w^/wo, 
the increase with temperature is faster since less energy 
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is required to excite the internal modes of the system and 
more states are available, correspondingly increasing S. 
Again, for larger values of ujr/uJo, the center of mass is 
the only active mode at low temperature and the division 
by N in the plots explains the lower value of S for larger 
N. Only at substantially larger temperatures will both 
internal and center of mass modes become active. The 
effect of symmetry is not very pronounced, although it is 
noticeable that the entropy is larger for fermionic than 
for bosonic systems due to the larger degeneracy. 

An interesting feature of the fermion plot is that for 
= 8 the limit of S/N for T — >■ is non- vanishing. 
The number of available states is finite reflecting that the 
ground state itself is degenerate. This is seen by simply 
counting of oscillator degeneracy for a two-dimensional 
harmonic oscillator. The lowest three quantum levels 
can hold 6 identical fermions (1,2, and 3, respectively), 
and the fourth can hold additionally four particles. This 
means that eight particles only occupy half of the last 
level, leaving the ground state as six times degenerate. 
This beavior at zero temperature then nicely indicate the 
presence of shell structure. 



C. Heat capacity and compressibility 

We now consider some derived thermodynamic quan- 
tities that are of experimental interest in many fields of 
physics. The heat capacity and the compressiblity are 
two such quantities that can be obtained from second 
derivatives of the partition function as discussed above. 
In Fig. ini we show the heat capacity per particle, Cy /N, 
defined in Eq.® as function of temperature for both 
bosons and fermions. They both start with an initial ac- 
tivation of the external trap mode, since the degenerate 
relative degrees of freedom all require higher temperature 
to be excited. After a delay, these internal modes are acti- 
vated, and the heat capacity increases with temperature. 
The delay and the rate of increase depend strongly on 
LUr/t-jQ with a slower increase of Cy for larger interaction 
ratios. 

The tendency to increase slower and in steps is re- 
lated to gaps in the energy spectrum. At high excita- 
tion energy, the spectrum becomes denser, and gap sizes 
larger than the temperature cannot appear. However, 
the presence of a gap in the low-energy spectrum is im- 
portant for the heat capacity in general. In fact, this is 
clearly seen by the fact that bosons have a flat profile 
for a region of low temperature, while the fermions have 
two flat plateaus. Fermions rise faster initially, and al- 
ways approach the equipartition heat capacity from be- 
low [Cv/N 2 for two dimensions). The larger gap 
causes a delay in the boson systems. The heat capacity 
slightly overshoots the equipartition value, oscillate back 
below the equipartion value (at a temperature outside 
the scale in Fig. [B]), and eventually approach the limit 
from below. The curve marked COM shows the heat ca- 
pacity when assuming that only the center of mass mode 
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FIG. 6: Heat capacities divided by NkB in two dimensions 
for boson (upper) and fermions (lower) systems, identified by 
the number of particles and frequency ratios, {N , u)r / too) ■ The 
black line shows the unsealed value for the mean field center 
of mass mode. 



is excited. All the other curves will approach this at very 
large temperatures when the internal structure is washed 
out. However, as the plots clearly show, the approach is 
very different for different particle numbers and interac- 
tion strengths. 

In Fig. [71 we show the isothermal compressibility per 
particle, kt/N, as defined in Ea. (fT5)) . The high temper- 
ature behavior of a harmonic system is 1/T. The com- 
pressiblity shows a small increase through a maximum at 
very small temperature followed by steady decrease to- 
wards zero at large temperature with a 1/T slope. The 
compressiblity indicates how easy it is to squeeze the 
system and a large value indicates that the system is 
very susceptible to compression. We see that the more 
strongly interacting systems (larger w^/wo) have larger 
kt- This comes from the fact that the attraction in these 
systems will make it energetically favourable to contract, 
due to the larger degeneracy of the interaction frequency, 
ujr- This also explains why kt increases with N. Com- 
paring bosons and fermions we find the same qualitative 
behavior. The fermions are slightly less compressible. 
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FIG. 7: Isothermal compressibility in two dimensions (I/Bt) 
divided by A'^ for boson (upper) and fermions (lower) systems, 
identified by the number of particles and frequency ratios, 

(A'', Ldr/tOo). 
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FIG. 8: Isothermal compressibility (I/Bt) divided by TV 
for fermionic systems in three dimensions, identified by the 
number of particles and frequency ratios, {N,u)r/u)o). 



with a sharper dependence on particle number and inter- 
action strength. This is a common feature of Fermi sys- 
tems and usually attributed to the Pauli principle. For 
comparison. Fig. |8] shows the compressibility for fermions 
in three spatial dimensions, which shows the same quali- 
tative behavior. The results for bosons are similar and we 
do not show them here. Notice that the peak features in 
the compressibility are sharper in the three-dimensional 
case. Our system is similar to an attractively interact- 
ing Fermi gas where superfluidity is expected to show 
a signature in the compressibility [T6j . Our results are 
consistent with the fact that phase transitions are less 
pronounced in lower dimensions. 



IV. SUMMARY AND CONCLUSIONS 

The harmonic approximation is extremely useful be- 
cause its simplicity allows transparent calculations of oth- 
erwise complicated properties. The only approximation 
lies in the choice of harmonic potentials acting on each 
particle and between pairs of particles. One can therefore 
think of the harmonic approximation as an effective inter- 
action scheme and subsequently investiagate the behavior 
of its predictions under changes in the input parameters. 
The latter should naturally be connected to whatever re- 
alistic physical system one is interested in studying. 

In this paper, we explore the harmonic approximation 
scheme by considering the Hamiltonian for a many-body 
system consisting of a given number of identical particles. 
In particular, we explore the consequences of symmetry 
requirements on the properties of the system. This is 
done with thermodynamic applications in mind since this 
is a venue where the quantum statistics plays a decisive 
role at low temperatures. We develop a new method 
for counting the number of correctly symmetrized many- 
body states that reduces the complications that arise from 
this fundamental problem of statistical mechanics and 
thermodynamics. 

The advantage of the harmonic approach is obviously 
that the energy spectrum is analytically known. How- 
ever, the degeneracy of each many-body state of given 
total energy still remains to be determined before the 
partition function is fully defined and possible to com- 
pute numerically. We design a novel procedure to obtain 
the degeneracy of each state, subject to requirements of 
symmetry and antisymmetry appropriate for bosons and 
fermions, respectively. We separate the completely per- 
mutation symmetric center of mass motion from the rela- 
tive motion, which then has to carry the symmetry corre- 
sponding to bosons or fermions. We count by subtracting 
the number of states of different quanta in the center of 
mass motion from the total number of non-interacting 
states of given symmetry. 

To demonstrate the method, we consider the case of a 
two-dimensional system with identical bosons or fermions 
(with no internal degrees of freedom) . Within the canon- 
ical ensemble, we compute the partition function, and 
from it the free energy, entropy, heat capacity, and com- 
pressibility of the system. This is done for a relatively 
small number of particles (up to 20). The method is a 
considerable improvement over the brute force method 
where one explicitly checks for symmetry properties by 
exchanging all pairs of particles one by one. However, it 
is still computationally involved when going beyond the 
particle numbers considered here. However, as is known 
from for instance the virial expansion [30|, it is often 
enough to consider small particle numbers and then ex- 
trapolate to large system sizes from this information. 

The effective harmonic interaction can be related to 
quantities in realistic systems and we discuss a case of 
great usefulness within the realm of ultracold atomic 
gases, that of two particles interaction in a harmonic 
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trap through a two-body interaction of zero-range. While 
bosons can interact in the s-wave channel originally con- 
sidered by Busch et al. [10], identical fermions must have 
an antisymmetric relative wave function and thus an s- 
wave interaction of short-range will vanish. We therefore 
have to consider the p-wave channel, but we find that the 
effective harmonic interaction frequencies are very similar 
to the s-wave case in the two-dimensional setup that we 
consider here. Therefore we have chosen to parametrize 
the discussion of the thermodynamic quantities by the 
harmonic oscillator frequency of the two-body interac- 
tion itself. One can then make the connection to a real- 
istic system by working backwards through the model of 
Busch et al.. 

Our numerical results show that the low-temperature 
behavior reflects shell structure for the particle numbers 
we study. The large-temperature equipartition limits are 
recovered for energy and heat capacity, and the transition 
from small to large temperature is fastest for smaller two- 
body interaction strength. The qualitative behavior is 
rather similar for bosons and fermions. However, bosons 
always have a gap in the low-energy spectrum since one 
quantum of excitation of the relative motion must be 
antisymmetric and therefore forbidden. This leads to a 
slower variation with temperature, since this gap has to 



be overcome before the number of available states goes 
up. Our results for the density of states turn out to scale 
with the number of relative excitation quanta in a manner 
that is very similar to the treatment of the many-body 
density of states for a uniform Fermi gas. Surprisingly, 
the bosonic many-body level density scales similarly to 
the fermionic one, although the particle number enters 
differently in our interpolation formulas. 

In future studies it will be interesting to consider 
also multi-component systems, something which is sim- 
ply done within the harmonic approach since the level 
counting can be factorized in the different components. 
Also, a study of the virial expansion based on the har- 
monic approach is currently on-going, both in two and 
three dimensions. In addition, we note that one dimen- 
sional system have attracted a lot of attention recently 
due to their realization in ultracold atomic physics Q . It 
would be interesting to test our prediction against some 
of the models that are being explored in the experiments 
and for which a number of exactly solvable many-body 
models are known Furthermore, recent experiments 
studying ultracold few-body two-component Fermi sys- 
tems (particle numbers of ten or less) Q would be an 
interesting comparison for the harmonic approximation. 
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